%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%``The Costs and Environmental Justice Concerns of...
%%... NIMBY in Solid Waste Disposal''*
%%% Phuong Ho - phuong.ho@snf.no %%%
%% Centre for Applied Research at Norwegian School of Economics %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This file has 3 parts.
%%% Part 1: Estimate the structural model using the main estimator
%%% Part 2: Calculate counterfactuals for universal NIMBY
%%% Part 3: Calculate counterfactuals for NIMBY in a single county
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
tic
cd 'C:\XYZ' %change working dir
cd matlab
diary(strcat('nimby_estimate_log_',date))
%%%% Import data 
mysample=readtable(strcat('input\nimby_estimate.csv'));
En=mysample(:,'realprice'); En=table2array(En);
Ex=mysample(:,{'dridis','localf'}); Ex=table2array(Ex); KEx=size(Ex,2);
Ins=mysample(:,{'sumothmktsizes50s','isumothmktsizes50s'});
Ins=table2array(Ins); KIns=size(Ins,2);
T=mysample(:,'T');T=table2array(T);
Y=mysample(:,'year');Y=table2array(Y);
Q=mysample(:,'quarter');Q=table2array(Q);
TC=mysample(:, 'TC'); TC=table2array(TC); TCmax=max(TC);
C=mysample(:,'C'); C=table2array(C); 
CJ=mysample(:,'CJ');CJ=table2array(CJ);
J=mysample(:,'J'); J=table2array(J);
F=mysample(:,'F'); F=table2array(F);
D=mysample(:,'D'); D=table2array(D);
TD=mysample(:,'TD'); TD=table2array(TD);
CD=mysample(:,'CD'); CD=table2array(CD);
TJ=mysample(:,'TJ'); TJ=table2array(TJ);
TF=mysample(:,'TF'); TF=table2array(TF);
q=mysample(:,'q'); q=table2array(q); 
local=mysample(:,'localf'); local=table2array(local);
MM=mysample(:,'mktsize'); MM=table2array(MM); %standard mktsize
M0=mysample(:,'cntyq0'); M0=table2array(M0); 
outq0=mysample(:,'outq0'); outq0=table2array(outq0);
q1=q>0; q0=q==0; lq=log(q); lq(q0)=0;
N=size(q,1); %number of observations
%%%% Get obs indices of each group in TC. Each cell is a tc group
uTC=unique(TC,'stable'); L=length(uTC); TC_idx=cell(L,1); vec=1:N; TC_rep=zeros(1,L);
TC_l=NaN(N,1);
for i=1:L
    %all obs indices of group tc=u(i). each cell is for 1 group
    TC_idx{i}=vec(TC==uTC(i)); 
    %representative indices of all groups, length(v)=length(u)
    %first representative index of each group tc=u(i)
    TC_rep(i)=TC_idx{i}(1,1); 
    TC_l(TC_idx{i})=size(TC_idx{i},2);
end
clear i vec l;
TCmj_idx=cell(N,1); %obs indices of k in TC_{j} but exclude facility j
for i=1:N, TCmj_idx{i}=setdiff(TC_idx{uTC==TC(i)},i); end
%%% Get a representative index in each group in J
u=unique(J,'stable'); L=length(u); J_idx=cell(L,1); vec=1:N; J_rep=zeros(1,L);
for i=1:L, J_idx{i}=vec(J==u(i)); J_rep(i)=J_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in T
u=unique(T,'stable'); L=length(u); T_idx=cell(L,1); vec=1:N; T_rep=zeros(1,L);
for i=1:L, T_idx{i}=vec(T==u(i)); T_rep(i)=T_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in C
u=unique(C,'stable'); L=length(u); C_idx=cell(L,1); vec=1:N; C_rep=zeros(1,L);
for i=1:L, C_idx{i}=vec(C==u(i)); C_rep(i)=C_idx{i}(1,1); end; clear i vec l u;
u=unique(CJ,'stable'); L=length(u); CJ_idx=cell(L,1); vec=1:N; CJ_rep=zeros(1,L);
for i=1:L, CJ_idx{i}=vec(CJ==u(i)); CJ_rep(i)=CJ_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TJ
u=unique(TJ,'stable'); L=length(u); TJ_idx=cell(L,1); vec=1:N; TJ_rep=zeros(1,L);
for i=1:L, TJ_idx{i}=vec(TJ==u(i)); TJ_rep(i)=TJ_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TF
u=unique(TF,'stable'); L=length(u); TF_idx=cell(L,1); vec=1:N; TF_rep=zeros(1,L);
for i=1:L, TF_idx{i}=vec(TF==u(i)); TF_rep(i)=TF_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in F
u=unique(F,'stable'); L=length(u); F_idx=cell(L,1); vec=1:N; F_rep=zeros(1,L);
for i=1:L, F_idx{i}=vec(F==u(i)); F_rep(i)=F_idx{i}(1,1); end; clear i vec l u;

M=M0; sj_q=my_accumarray(TC_idx,q); %sum of q(j) within TC (market q)
sh=q./(M); outq=outq0; outsh=outq./(M); %outside option's share
oq1=outq>0; oq0=outq==0;
NN=length(q);
sqw=q./NN./mean(M); outqw=outq./NN./mean(M); outqwrep=outqw(TC_rep);
outshrep=outsh(TC_rep);
%%%% Building covariates
%Declare type of fixed effects
XFE=[F T];
%Create dummies for FEs.
XFD=dummy_full(XFE);
XFD(:,size(XFD,2))=[]; %remove the last column as the base

%%%%%%%%%%%%%%% PART 1: ESTIMATE THE FIRST ESTIMATOR %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% SPECIFICATION 1: J FE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X0=[En Ex]; KX0=size(X0,2);
X0scl=zeros(size(X0)); X0scl=X0 - repmat(mean(X0,1), N,1);
sclX=std(X0,1); X0scl=X0scl ./ repmat(sclX,N,1);
X=[X0scl XFD]; KX=size(X,2); %true X covariates for computation
%%%% NOW, ESTIMATE THE STRUCTURAL MODEL %%%%%%%
beta0=zeros(KX,1);
fun_beta=@(beta)LLHs(beta,X,TC_idx,q0,q1,sqw);
disp('Specification 1: Facility fixed effects. List of Iterations');
ops=optimoptions('fminunc','Algorithm','trust-region','Display','iter',...
    'MaxIterations',400,'SpecifyObjectiveGradient',true,...
    'StepTolerance',1e-7,'FunctionTolerance',1e-7,'OptimalityTolerance',1e-6); 
[beta,feval,exitflag,output,grad,hessian]=fminunc(fun_beta,beta0,ops);
disp('Specification 1: Facility fixed effects. Point estimates are');
beta(1:KX0)=beta(1:KX0)./sclX';%descaled beta
fprintf('Estimate of the first estimator (controlling fixed effects)');
beta(1:KX0)

%%% Calculate predicted shares
X=[En Ex XFD]; KX=size(X,2); %descaled X
betaX=beta;
[pr,CS]=pr_fun(betaX,X,TC_idx); %counterfactual market share
el=beta(1).*En.*(1-pr); %price elasticity
eldis=beta(2).*Ex.*(1-pr); %distance elasticity
disp('price elasticity');
el'*q./sum(q) % price elasticity
disp('transport elasticity');
eldis'*q./sum(q) % transport elasticity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% calculate status quo marginal cost %%%%%%%%%%%%%%%%%%%%%%%%
dsdp=beta(1)*pr.*(1-pr); %dsj/dpj
nume=M./mean(M).*pr; 
nume=my_accumarray(TF_idx,nume); %sum of shares over county for a given facility
deno=my_accumarray(TF_idx,M./mean(M).*dsdp);
%assume a given facility sets same price for all counties
%then we only sum shares and dsdp over counties for a facility
pmc=-nume./deno; 
mc=En-pmc;
markup=pmc./(En-pmc); % this is (p-mc)/mc
disp('average markup');
markup'*q./sum(q)  % markup
qj=my_accumarray(TF_idx,M./100000.*pr); %total trash going to facility j at t
%reduce dimension CTJx1 to TJx1. Reduce dimension using TF_rep
qj_jt=qj(TF_rep); lnQj_jt=log(qj_jt); lnQj_jt2=log(qj_jt).^2;
qj2_jt=qj_jt.^2; qj3_jt=qj_jt.^3;
mc_jt=mc(TF_rep); %can test by check: [TJ(1:30) mc(1:30)]
disp('cost coefficient: mc=gamma_1*Q + gamma_2*Q^2');
F_jt=F(TF_rep); T_jt=T(TF_rep);
regression2(mc_jt,[qj_jt qj2_jt],[],[],[F_jt T_jt]); clear Bdel;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% SPECIFICATION 2: IV + linear control function + J FE %%%%%%%
[Bdel,R]=regression2(En,[Ins Ex],[],[],[F T]); clear Bdel;
X0=[En Ex R]; KX0=size(X0,2);
X0scl=zeros(size(X0)); X0scl=X0 - repmat(mean(X0,1), N,1);
sclX=std(X0,1); X0scl=X0scl ./ repmat(sclX,N,1);
X=[X0scl XFD]; KX=size(X,2); %true X covariates for computation
beta0=zeros(KX,1);
fun_beta=@(beta)LLHs(beta,X,TC_idx,q0,q1,sqw);
disp('Specification 2: linear control + facility FE. List of iterations:');
ops=optimoptions('fminunc','Algorithm','trust-region','Display','iter',...
    'MaxIterations',400,'SpecifyObjectiveGradient',true,...
    'StepTolerance',1e-7,'FunctionTolerance',1e-7,'OptimalityTolerance',1e-6); 
[beta,feval,exitflag,output,grad,hessian]=fminunc(fun_beta,beta0,ops);
disp('Specification 2: linear control + facility FE. Estimates are');
beta(1:KX0)=beta(1:KX0)./sclX';%descaled beta
fprintf('Second estimator (linear control function with fixed effects)');
beta(1:KX0)

%%% Calculate predicted shares
X=[En Ex R XFD]; KX=size(X,2); %descaled X
betaX=beta;
[pr,CS]=pr_fun(betaX,X,TC_idx); %counterfactual market share
el=beta(1).*En.*(1-pr); %price elasticity
eldis=beta(2).*Ex.*(1-pr); %distance elasticity
disp('price elasticity');
el'*q./sum(q) % price elasticity
disp('transport elasticity');
eldis'*q./sum(q) % transport elasticity
save(strcat('output\result_main_estimator'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% calculate status quo marginal cost %%%%%%%%%%%%%%%%%%%%%%%%
dsdp=beta(1)*pr.*(1-pr); %dsj/dpj
nume=M./mean(M).*pr;
nume=my_accumarray(TF_idx,nume); 
deno=my_accumarray(TF_idx,M./mean(M).*dsdp);
pmc=-nume./deno; 
mc=En-pmc;
markup=pmc./(En-pmc); % this is (p-mc)/mc
disp('average markup');
markup'*q./sum(q)  % markup
qj=my_accumarray(TF_idx,M./100000.*pr); %total trash going to facility j
%reduce dimension CTJx1 to TJx1. Reduce dimension using TF_rep
qj_jt=qj(TF_rep); lnQj_jt=log(qj_jt); lnQj_jt2=log(qj_jt).^2;
qj2_jt=qj_jt.^2; qj3_jt=qj_jt.^3;
mc_jt=mc(TF_rep); %can test by check: [TJ(1:30) mc(1:30)]
F_jt=F(TF_rep); T_jt=T(TF_rep);
regression2(mc_jt,[qj_jt qj2_jt],[],[],[F_jt T_jt]); clear Bdel;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% SPECIFICATION 3: IV + quadratic control function + J FE %%%%%%%
R2=R.^2;
X0=[En Ex R R2]; KX0=size(X0,2);
X0scl=zeros(size(X0)); X0scl=X0 - repmat(mean(X0,1), N,1);
sclX=std(X0,1); X0scl=X0scl ./ repmat(sclX,N,1);
X=[X0scl XFD]; KX=size(X,2); %true X covariates for computation
beta0=zeros(KX,1);
fun_beta=@(beta)LLHs(beta,X,TC_idx,q0,q1,sqw);
disp('Specification 3: quadratic control + facility FE. List of iterations:');
ops=optimoptions('fminunc','Algorithm','trust-region','Display','iter',...
    'MaxIterations',400,'SpecifyObjectiveGradient',true,...
    'StepTolerance',1e-7,'FunctionTolerance',1e-7,'OptimalityTolerance',1e-6); 
[beta,feval,exitflag,output,grad,hessian]=fminunc(fun_beta,beta0,ops);
disp('Specification 2: linear control + facility FE. Estimates are');
beta(1:KX0)=beta(1:KX0)./sclX';%descaled beta
fprintf('Third estimator (quadratic control function with fixed effects)');
beta(1:KX0)
beta(KX0)*10^6

%%% Calculate predicted shares
X=[En Ex R R2 XFD]; KX=size(X,2); %descaled X
betaX=beta;
[pr,CS]=pr_fun(betaX,X,TC_idx); %counterfactual market share
el=beta(1).*En.*(1-pr); %price elasticity
eldis=beta(2).*Ex.*(1-pr); %distance elasticity
disp('price elasticity');
el'*q./sum(q) % price elasticity
disp('transport elasticity');
eldis'*q./sum(q) % transport elasticity

dsdp=beta(1)*pr.*(1-pr); %dsj/dpj
nume=M./mean(M).*pr;
nume=my_accumarray(TF_idx,nume);
deno=my_accumarray(TF_idx,M./mean(M).*dsdp); 
pmc=-nume./deno; 
mc=En-pmc;
markup=pmc./(En-pmc); % this is (p-mc)/mc
markup'*q./sum(q)  % markup
qj=my_accumarray(TF_idx,M./100000.*pr); %total trash going to facility j
qj_jt=qj(TF_rep); lnQj_jt=log(qj_jt); lnQj_jt2=log(qj_jt).^2;
qj2_jt=qj_jt.^2; qj3_jt=qj_jt.^3;
mc_jt=mc(TF_rep); 
F_jt=F(TF_rep); T_jt=T(TF_rep);
regression2(mc_jt,[qj_jt qj2_jt],[],[],[F_jt T_jt]); clear Bdel;

%%%%%%%%%%%% ENDING: ESTIMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%
%%%%%%%%% PART 3: PREDICTION AND COUNTERFACTUAL CALCULATION %%%%%%%%
%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% SAVE FITTED VALUES %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%load the main estimates
load(strcat('output\result_main_estimator'));
M=sj_q; dsdp=beta(1)*pr.*(1-pr); %dsj/dpj
nume=M./mean(M).*pr;
nume=my_accumarray(TF_idx,nume); 
deno=my_accumarray(TF_idx,M./mean(M).*dsdp); 
pmc=-nume./deno; 
mc=En-pmc;
markup=pmc./(En-pmc);
qj=my_accumarray(TF_idx,M./100000.*pr); 
qj_jt=qj(TF_rep); lnQj_jt=log(qj_jt);
qj2_jt=qj_jt.^2; qj3_jt=qj_jt.^3;
mc_jt=mc(TF_rep); 
[c,bint,c_jt]=regress(mc_jt,[qj_jt qj2_jt]);
c_cjt=zeros(N,1); c_cjt(TF_rep)=c_jt; %input value
c_cjt=my_accumarray(TF_idx,c_cjt); %spread values, (since other values are 0)
Bq=pr.*M; %fitted status quo trash amount in a facility
mysample.sj_q=sj_q; %use sj_q because this is the symbol in STATA
mysample.pr=pr; %fitted share
mysample.Bq=Bq;
mysample.el=el;
mysample.CS=CS;
mysample.pmc=pmc;
mysample.markup=markup;
save(strcat('output\result_main_estimator'));
writetable(mysample,'output\nimby_statusquo.xlsx','Sheet',1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% COUNTERFACTUAL WITHOUT PRICE ADJUSTMENT %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% IMPORT BAN. WITHOUT PRICE ADJUSTMENT %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(strcat('output\result_main_estimator'));
M=sj_q;
 % import new data, calculate new q, CS given estimated beta and new data
TC(mysample.localf==0)=99999;
%% Get obs indices of each group in TC. Each cell is a tc group
uTC=unique(TC); L=length(uTC); TC_idx=cell(L,1); vec=1:N; TC_rep=zeros(1,L);
TC_l=NaN(N,1);
for i=1:L
    TC_idx{i}=vec(TC==uTC(i)); 
    TC_rep(i)=TC_idx{i}(1,1); 
    TC_l(TC_idx{i})=size(TC_idx{i},2);
end
 %%calculate counterfactual share and q and policy variables
X=[En Ex R XFD]; KX=size(X,2); %descaled X
[Cpr,CCS]=pr_fun(beta,X,TC_idx); %counterfactual market share
Cq=Cpr.*M; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=En;
mysample(mysample.localf==0,:)=[];
writetable(mysample,'output\nimby_counterfactual_importBan.xlsx','Sheet',1);

%%%%%%%%%%%%% IMPORT TAX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(strcat('output\result_main_estimator'));
N=size(X,1); %number of observations
NEn=En; Cpr=zeros(N,1); Cq=zeros(N,1); CCS=zeros(N,1);
NEn(mysample.localf==0)=En(mysample.localf==0).*1.10; %regulated price
CX=[NEn Ex R XFD];
[Cpr,CCS]=pr_fun(beta,CX,TC_idx); %counterfactual market share
 Cq=Cpr.*sj_q; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=En;%Cprice is new price facility set, excl tax
writetable(mysample,'output\nimby_counterfactual_importTax.xlsx','Sheet',1);

%%%%%%COUNTERFACTUAL: WITH PRICE ADJUSTMENT %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% IMPORT BAN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(strcat('output\result_main_estimator'));
M=sj_q; mM=mean(M);
%%%%%%%%%%% Counterfactuals %%%%%%%%%%%%%%%%%%%%%%%%%%%
%Claim variables and define counterfactual for fuel tax and import ban;
%For trash tax and import tax, modify the function that calculate eq price
%keep local observations for mysample. Define obs to remove
%i.e. remove nonlocalf and not 2010
foreign=mysample.localf==0 | mysample.year~=2010;% zeros(N,1);
mysample(foreign==1,:)=[]; %only keep local flow sample
c_cjt(foreign==1,:)=[];
M(foreign==1,:)=[];
NXFD=XFD(foreign==0,:); NR=R(foreign==0,:); 
%redefine variables from mysample
En=mysample(:,'realprice'); En=table2array(En);
Ex=mysample(:,{'dridis','localf'}); Ex=table2array(Ex); KEx=size(Ex,2);
Ins=mysample(:,{'sumothmktsizes50s','isumothmktsizes50s'});
Ins=table2array(Ins); KIns=size(Ins,2);
T=mysample(:,'T');T=table2array(T);
Y=mysample(:,'year');Y=table2array(Y);
Q=mysample(:,'quarter');Q=table2array(Q);
TC=mysample(:, 'TC'); TC=table2array(TC); TCmax=max(TC);
C=mysample(:,'C'); C=table2array(C); 
CJ=mysample(:,'CJ');CJ=table2array(CJ);
J=mysample(:,'J'); J=table2array(J);
F=mysample(:,'F'); F=table2array(F);
D=mysample(:,'D'); D=table2array(D);
TD=mysample(:,'TD'); TD=table2array(TD);
CD=mysample(:,'CD'); CD=table2array(CD);
TJ=mysample(:,'TJ'); TJ=table2array(TJ);
TF=mysample(:,'TF'); TF=table2array(TF);
q=mysample(:,'q'); q=table2array(q); 
q1=q>0; q0=q==0; lq=log(q); lq(q0)=0;
N=size(q,1); %number of observations
%%%% Get obs indices of each group in TC. Each cell is a tc group
uTC=unique(TC,'stable'); L=length(uTC); TC_idx=cell(L,1); vec=1:N; TC_rep=zeros(1,L);
TC_l=NaN(N,1);
for i=1:L
    %all obs indices of group tc=u(i). each cell is for 1 group
    TC_idx{i}=vec(TC==uTC(i)); 
    %representative indices of all groups, length(v)=length(u)
    %first representative index of each group tc=u(i)
    TC_rep(i)=TC_idx{i}(1,1); 
    TC_l(TC_idx{i})=size(TC_idx{i},2);
end
clear i vec l;
TCmj_idx=cell(N,1); %obs indices of k in TC_{j} but exclude facility j
for i=1:N, TCmj_idx{i}=setdiff(TC_idx{uTC==TC(i)},i); end
%%% Get a representative index in each group in J
u=unique(J,'stable'); L=length(u); J_idx=cell(L,1); vec=1:N; J_rep=zeros(1,L);
for i=1:L, J_idx{i}=vec(J==u(i)); J_rep(i)=J_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in T
u=unique(T,'stable'); L=length(u); T_idx=cell(L,1); vec=1:N; T_rep=zeros(1,L);
for i=1:L, T_idx{i}=vec(T==u(i)); T_rep(i)=T_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in C
u=unique(C,'stable'); L=length(u); C_idx=cell(L,1); vec=1:N; C_rep=zeros(1,L);
for i=1:L, C_idx{i}=vec(C==u(i)); C_rep(i)=C_idx{i}(1,1); end; clear i vec l u;
u=unique(CJ,'stable'); L=length(u); CJ_idx=cell(L,1); vec=1:N; CJ_rep=zeros(1,L);
for i=1:L, CJ_idx{i}=vec(CJ==u(i)); CJ_rep(i)=CJ_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TJ
u=unique(TJ,'stable'); L=length(u); TJ_idx=cell(L,1); vec=1:N; TJ_rep=zeros(1,L);
for i=1:L, TJ_idx{i}=vec(TJ==u(i)); TJ_rep(i)=TJ_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TF
u=unique(TF,'stable'); L=length(u); TF_idx=cell(L,1); vec=1:N; TF_rep=zeros(1,L);
for i=1:L, TF_idx{i}=vec(TF==u(i)); TF_rep(i)=TF_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in F
u=unique(F,'stable'); L=length(u); F_idx=cell(L,1); vec=1:N; F_rep=zeros(1,L);
for i=1:L, F_idx{i}=vec(F==u(i)); F_rep(i)=F_idx{i}(1,1); end; clear i vec l u;

%counterfactual values
CEn=zeros(N,1); Cpr=zeros(N,1); Cq=zeros(N,1); CCS=zeros(N,1);
Cpmc=zeros(N,1); Cmarkup=zeros(N,1); NEn=zeros(N,1);
CX=[En Ex NR NXFD];
%%%%%%%%%%%% Restrict the data to year 2010 %%%%%%%%%%%%%%%%%%
%%Run estimation for every quarter in the year %%%
for iQ=1:4
%iQ=1;
S=find(Y==2010&Q==iQ); %index for counterfactual parts
%%% Only care about the experiment period S.
%Vectors starting with S have dimension Sx1 %%
Sc_cjt=c_cjt(S); SM=M(S);
STC=TC(S); STF=TF(S); SEn=En(S); SN=size(STC,1);
SX=CX(S,:); KSX=size(SX,2);
beta1=betaX(1); KbetaX=size(betaX,1);
eX2beta2=exp(SX(:,2:KSX)*betaX(2:KbetaX)); 
%%% Get a representative index in each group in TJ
u=unique(STF,'stable'); L=length(u); STF_idx=cell(L,1); vec=1:SN; STF_rep=zeros(1,L);
for i=1:L, STF_idx{i}=vec(STF==u(i)); STF_rep(i)=STF_idx{i}(1,1); end;
clear i vec l u;
%%% Get a representative index in each group in TJ
u=unique(STC,'stable'); L=length(u); STC_idx=cell(L,1); vec=1:SN; STC_rep=zeros(1,L);
for i=1:L, STC_idx{i}=vec(STC==u(i)); STC_rep(i)=STC_idx{i}(1,1); end; clear i vec l u;
%%Shorten price and cost to TJx1, instead of CTJx1
P_jt=SEn(STF_rep); Sc_jt=Sc_cjt(STF_rep);
%old price/observed price is guessing point
%%Define the equation to solve for equilibrium price
fun_price_importban=@(P_jt)solveprice_importban(beta1,P_jt,eX2beta2,Sc_jt,c,SM,STC_idx,STF_idx,STF_rep,SN,mM);
opfsolve=optimoptions('fsolve','Algorithm','levenberg-marquardt','Display','iter',...
    'MaxIterations',400,'MaxFunctionEvaluations',1e15,...
    'StepTolerance',1e-10,'FunctionTolerance',1e-10,'OptimalityTolerance',1e-6); 
%%Calculate eq price for trash tax
[SP,feval,exitflag,output,jacob]=fsolve(fun_price_importban,P_jt,opfsolve);
%span P_jt to P_cjt
CEn(intersect(TF_rep,S))=SP; %P_jt in period S
CEn(S)=my_accumarray(STF_idx,CEn(S));%span. same as egen x=mean(.) or total
SX=[CEn(S) SX(:,2:KSX)];
[Cpr(S),CCS(S)]=pr_fun(betaX,SX,STC_idx); %counterfactual market share
%% calculate new price cost margin
dsdp(S)=betaX(1)*Cpr(S).*(1-Cpr(S)); %dsj/dpj
nume(S)=SM./mean(M).*Cpr(S);
nume(S)=my_accumarray(STF_idx,nume(S));
deno(S)=my_accumarray(STF_idx,nume(S).*dsdp(S));
Cpmc(S)=-nume(S)./deno(S); Cmarkup(S)=Cpmc(S)./(CEn(S)-Cpmc(S));
end
%out of the loop S
Cq=Cpr.*M; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=CEn; %new eq price, excl tax. Nrealprice=CEn*tax...
mysample.Cpmc=Cpmc; mysample.Cmarkup=Cmarkup;
writetable(mysample,'output\nimby_counterfactual_importBan_with_eqprice.xlsx','Sheet',1);

%%%%%%%%%%%%% IMPORT TAX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(strcat('output\result_main_estimator'));
M=sj_q; mM=mean(M);
%%%%%%%%%%% Counterfactuals %%%%%%%%%%%%%%%%%%%%%%%%%%%
%Claim variables and define counterfactual for fuel tax and import ban;
%For trash tax and import tax, modify the function that calculate eq price
%keep local observations for mysample. Define obs to remove
%i.e. remove nonlocalf and not 2010
remove=mysample.year~=2010;% zeros(N,1);
mysample(remove==1,:)=[]; %only keep local flow sample
c_cjt(remove==1,:)=[];
M(remove==1,:)=[];
NXFD=XFD(remove==0,:); NR=R(remove==0,:); 
%redefine variables from mysample
En=mysample(:,'realprice'); En=table2array(En);
Ex=mysample(:,{'dridis','localf'}); Ex=table2array(Ex); KEx=size(Ex,2);
T=mysample(:,'T');T=table2array(T);
Y=mysample(:,'year');Y=table2array(Y);
Q=mysample(:,'quarter');Q=table2array(Q);
TC=mysample(:, 'TC'); TC=table2array(TC); TCmax=max(TC);
C=mysample(:,'C'); C=table2array(C); 
F=mysample(:,'F'); F=table2array(F);
TF=mysample(:,'TF'); TF=table2array(TF);
q=mysample(:,'q'); q=table2array(q); 
local=mysample(:,'localf'); local=table2array(local);
q1=q>0; q0=q==0; lq=log(q); lq(q0)=0;
N=size(q,1); %number of observations
%%%% Get obs indices of each group in TC. Each cell is a tc group
uTC=unique(TC,'stable'); L=length(uTC); TC_idx=cell(L,1); vec=1:N; TC_rep=zeros(1,L);
TC_l=NaN(N,1);
for i=1:L
    %all obs indices of group tc=u(i). each cell is for 1 group
    TC_idx{i}=vec(TC==uTC(i)); 
    %representative indices of all groups, length(v)=length(u)
    %first representative index of each group tc=u(i)
    TC_rep(i)=TC_idx{i}(1,1); 
    TC_l(TC_idx{i})=size(TC_idx{i},2);
end
clear i vec l;
TCmj_idx=cell(N,1); %obs indices of k in TC_{j} but exclude facility j
for i=1:N, TCmj_idx{i}=setdiff(TC_idx{uTC==TC(i)},i); end
%%% Get a representative index in each group in T
u=unique(T,'stable'); L=length(u); T_idx=cell(L,1); vec=1:N; T_rep=zeros(1,L);
for i=1:L, T_idx{i}=vec(T==u(i)); T_rep(i)=T_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in C
u=unique(C,'stable'); L=length(u); C_idx=cell(L,1); vec=1:N; C_rep=zeros(1,L);
for i=1:L, C_idx{i}=vec(C==u(i)); C_rep(i)=C_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TF
u=unique(TF,'stable'); L=length(u); TF_idx=cell(L,1); vec=1:N; TF_rep=zeros(1,L);
for i=1:L, TF_idx{i}=vec(TF==u(i)); TF_rep(i)=TF_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in F
u=unique(F,'stable'); L=length(u); F_idx=cell(L,1); vec=1:N; F_rep=zeros(1,L);
for i=1:L, F_idx{i}=vec(F==u(i)); F_rep(i)=F_idx{i}(1,1); end; clear i vec l u;

%%%%%%%%%%%%%%%%%%%%%%%%below is fixed %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CEn=zeros(N,1); Cpr=zeros(N,1); Cq=zeros(N,1); CCS=zeros(N,1);
Cpmc=zeros(N,1); Cmarkup=zeros(N,1); NEn=zeros(N,1);
CX=[En Ex NR NXFD];
%%%%%%%%%%%% Restrict the data to year 2010 %%%%%%%%%%%%%%%%%%
%%Run estimation for every quarter in the year %%%
for iQ=1:4
S=find(Y==2010&Q==iQ); %define counterfactual period
%%% Only care about the experiment period S.
%Vectors starting with S have dimension Sx1 %%
Sc_cjt=c_cjt(S); SM=M(S)./1; 
STC=TC(S); STF=TF(S); SEn=En(S); SN=size(STC,1);
SX=CX(S,:); KSX=size(SX,2); Slocal=local(S,:);
beta1=betaX(1); KbetaX=size(betaX,1);
eX2beta2=exp(SX(:,2:KSX)*betaX(2:KbetaX));
%%% Get a representative index in each group in TJ
u=unique(STF,'stable'); L=length(u); STF_idx=cell(L,1); vec=1:SN; STF_rep=zeros(1,L);
for i=1:L, STF_idx{i}=vec(STF==u(i)); STF_rep(i)=STF_idx{i}(1,1); end;
clear i vec l u;
%%% Get a representative index in each group in TJ
u=unique(STC,'stable'); L=length(u); STC_idx=cell(L,1); vec=1:SN; STC_rep=zeros(1,L);
for i=1:L, STC_idx{i}=vec(STC==u(i)); STC_rep(i)=STC_idx{i}(1,1); end; clear i vec l u;
%%Shorten price and cost to TJx1, instead of CTJx1
P_jt=SEn(STF_rep); Sc_jt=Sc_cjt(STF_rep);
%%Define the equation to solve for equilibrium price
fun_price_importtax=@(P_jt)solveprice_importtax(beta1,P_jt,eX2beta2,Sc_jt,c,SM,Slocal,STC_idx,STF_idx,STF_rep,SN,mM);
opfsolve=optimoptions('fsolve','Algorithm','levenberg-marquardt','Display','iter',...
    'MaxIterations',400,'MaxFunctionEvaluations',1e15,...
    'StepTolerance',1e-12,'FunctionTolerance',1e-12,'OptimalityTolerance',1e-5); 
%%Calculate eq price for trash tax
[SP,feval,exitflag,output,jacob]=fsolve(fun_price_importtax,P_jt,opfsolve);
%span P_jt to P_cjt
CEn(intersect(TF_rep,S))=SP; %P_jt;
CEn(S)=my_accumarray(STF_idx,CEn(S));
NEn(S)=CEn(S);
NEn(S & local(S)==0)=CEn(S & local(S)==0)*1.10; % price includes tax
SX=[NEn(S) SX(:,2:KSX)];
[Cpr(S),CCS(S)]=pr_fun(betaX,SX,STC_idx); %counterfactual market share
%% calculate new price cost margin
dsdp(S)=betaX(1)*Cpr(S).*(1-Cpr(S)); %dsj/dpj
nume(S)=SM./mean(M).*Cpr(S);
nume(S)=my_accumarray(STF_idx,nume(S));
deno(S)=my_accumarray(STF_idx,nume(S).*dsdp(S));
Cpmc(S)=-nume(S)./deno(S); Cmarkup(S)=Cpmc(S)./(CEn(S)-Cpmc(S));
end
%out of the loop S
Cq=Cpr.*M; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=CEn; 
mysample.Cpmc=Cpmc; mysample.Cmarkup=Cmarkup;
writetable(mysample,'output\nimby_counterfactual_importTax_with_eqprice.xlsx','Sheet',1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% COUNTERFACTUAL: SINGLE COUNTY: MARIN %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% IMPORT BAN. WITHOUT PRICE ADJUSTMENT %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
counter=["Marin","Solano"];
for coun=1:2
load(strcat('output\result_main_estimator'));
M=sj_q;
%%%%%%%%%%% Counterfactuals %%%%%%%%%%%%%%%%%%%%%%%%%%%
%Claim variables and define counterfactual for fuel tax and import ban;
%For trash tax and import tax, modify the function that calculate eq price
%keep local observations for mysample. Define obs to remove
%i.e. remove nonlocalf and TO marin (import ban, not export ban) or all not 2010
remove=(mysample.localf==0&mysample.desco==counter(coun)) | mysample.year~=2010;
mysample(remove==1,:)=[]; %only keep local flow sample
c_cjt(remove==1,:)=[];
M(remove==1,:)=[];
XFD(remove==1,:)=[]; R(remove==1,:)=[]; 
%redefine variables from mysample
En=mysample(:,'realprice'); En=table2array(En);
Ex=mysample(:,{'dridis','localf'}); Ex=table2array(Ex); KEx=size(Ex,2);
T=mysample(:,'T');T=table2array(T);
Y=mysample(:,'year');Y=table2array(Y);
Q=mysample(:,'quarter');Q=table2array(Q);
TC=mysample(:, 'TC'); TC=table2array(TC); TCmax=max(TC);
C=mysample(:,'C'); C=table2array(C); 
CJ=mysample(:,'CJ');CJ=table2array(CJ);
J=mysample(:,'J'); J=table2array(J);
F=mysample(:,'F'); F=table2array(F);
TJ=mysample(:,'TJ'); TJ=table2array(TJ);
TF=mysample(:,'TF'); TF=table2array(TF);
q=mysample(:,'q'); q=table2array(q); 
local=mysample(:,'localf'); local=table2array(local);
q1=q>0; q0=q==0; lq=log(q); lq(q0)=0;
N=size(q,1); %number of observations
%%%% Get obs indices of each group in TC. Each cell is a tc group
uTC=unique(TC,'stable'); L=length(uTC); TC_idx=cell(L,1); vec=1:N; TC_rep=zeros(1,L);
TC_l=NaN(N,1);
for i=1:L
    %all obs indices of group tc=u(i). each cell is for 1 group
    TC_idx{i}=vec(TC==uTC(i)); 
    %representative indices of all groups, length(v)=length(u)
    %first representative index of each group tc=u(i)
    TC_rep(i)=TC_idx{i}(1,1); 
    TC_l(TC_idx{i})=size(TC_idx{i},2);
end
clear i vec l;
%%calculate counterfactual share and q and policy variables
X=[En Ex R XFD]; KX=size(X,2); %descaled X
[Cpr,CCS]=pr_fun(beta,X,TC_idx); %counterfactual market share
Cq=Cpr.*M; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=En;
writetable(mysample,strcat('output\nimby_counterfactual_importBan_',counter(coun),'.xlsx'),'Sheet',1);

%%%%%%%%%%%%% IMPORT TAX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(strcat('output\result_main_estimator'));
N=size(X,1); %number of observations
NEn=En; Cpr=zeros(N,1); Cq=zeros(N,1); CCS=zeros(N,1);
NEn(mysample.localf==0&mysample.desco==counter(coun))=En(mysample.localf==0&mysample.desco==counter(coun)).*1.10; %regulated price
CX=[NEn Ex R XFD];
[Cpr,CCS]=pr_fun(beta,CX,TC_idx); %counterfactual market share
 Cq=Cpr.*sj_q; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=En;
writetable(mysample,strcat('output\nimby_counterfactual_importTax_',counter(coun),'.xlsx'),'Sheet',1);

%%%%%%COUNTERFACTUAL: WITH PRICE ADJUSTMENT %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% IMPORT BAN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(strcat('output\result_main_estimator'));
M=sj_q; mM=mean(M);
%%%%%%%%%%% Counterfactuals %%%%%%%%%%%%%%%%%%%%%%%%%%%
%Claim variables and define counterfactual for fuel tax and import ban;
%For trash tax and import tax, modify the function that calculate eq price
%keep local observations for mysample. Define obs to remove
%i.e. remove nonlocalf and not 2010
remove=(mysample.localf==0&mysample.desco==counter(coun)) | mysample.year~=2010;% zeros(N,1);
mysample(remove==1,:)=[]; %only keep local flow sample
c_cjt(remove==1,:)=[];
M(remove==1,:)=[];
NXFD=XFD(remove==0,:); NR=R(remove==0,:); 
%redefine variables from mysample
En=mysample(:,'realprice'); En=table2array(En);
Ex=mysample(:,{'dridis','localf'}); Ex=table2array(Ex); KEx=size(Ex,2);
Ins=mysample(:,{'sumothmktsizes50s','isumothmktsizes50s'});
Ins=table2array(Ins); KIns=size(Ins,2);
T=mysample(:,'T');T=table2array(T);
Y=mysample(:,'year');Y=table2array(Y);
Q=mysample(:,'quarter');Q=table2array(Q);
TC=mysample(:, 'TC'); TC=table2array(TC); TCmax=max(TC);
C=mysample(:,'C'); C=table2array(C); 
CJ=mysample(:,'CJ');CJ=table2array(CJ);
J=mysample(:,'J'); J=table2array(J);
F=mysample(:,'F'); F=table2array(F);
D=mysample(:,'D'); D=table2array(D);
TD=mysample(:,'TD'); TD=table2array(TD);
CD=mysample(:,'CD'); CD=table2array(CD);
TJ=mysample(:,'TJ'); TJ=table2array(TJ);
TF=mysample(:,'TF'); TF=table2array(TF);
q=mysample(:,'q'); q=table2array(q); 
local=mysample(:,'localf'); local=table2array(local);
q1=q>0; q0=q==0; lq=log(q); lq(q0)=0;
N=size(q,1); %number of observations
%%%% Get obs indices of each group in TC. Each cell is a tc group
uTC=unique(TC,'stable'); L=length(uTC); TC_idx=cell(L,1); vec=1:N; TC_rep=zeros(1,L);
TC_l=NaN(N,1);
for i=1:L
    %all obs indices of group tc=u(i). each cell is for 1 group
    TC_idx{i}=vec(TC==uTC(i)); 
    %representative indices of all groups, length(v)=length(u)
    %first representative index of each group tc=u(i)
    TC_rep(i)=TC_idx{i}(1,1); 
    TC_l(TC_idx{i})=size(TC_idx{i},2);
end
clear i vec l;
TCmj_idx=cell(N,1); %obs indices of k in TC_{j} but exclude facility j
for i=1:N, TCmj_idx{i}=setdiff(TC_idx{uTC==TC(i)},i); end
%%% Get a representative index in each group in J
u=unique(J,'stable'); L=length(u); J_idx=cell(L,1); vec=1:N; J_rep=zeros(1,L);
for i=1:L, J_idx{i}=vec(J==u(i)); J_rep(i)=J_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in T
u=unique(T,'stable'); L=length(u); T_idx=cell(L,1); vec=1:N; T_rep=zeros(1,L);
for i=1:L, T_idx{i}=vec(T==u(i)); T_rep(i)=T_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in C
u=unique(C,'stable'); L=length(u); C_idx=cell(L,1); vec=1:N; C_rep=zeros(1,L);
for i=1:L, C_idx{i}=vec(C==u(i)); C_rep(i)=C_idx{i}(1,1); end; clear i vec l u;
u=unique(CJ,'stable'); L=length(u); CJ_idx=cell(L,1); vec=1:N; CJ_rep=zeros(1,L);
for i=1:L, CJ_idx{i}=vec(CJ==u(i)); CJ_rep(i)=CJ_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TJ
u=unique(TJ,'stable'); L=length(u); TJ_idx=cell(L,1); vec=1:N; TJ_rep=zeros(1,L);
for i=1:L, TJ_idx{i}=vec(TJ==u(i)); TJ_rep(i)=TJ_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TF
u=unique(TF,'stable'); L=length(u); TF_idx=cell(L,1); vec=1:N; TF_rep=zeros(1,L);
for i=1:L, TF_idx{i}=vec(TF==u(i)); TF_rep(i)=TF_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in F
u=unique(F,'stable'); L=length(u); F_idx=cell(L,1); vec=1:N; F_rep=zeros(1,L);
for i=1:L, F_idx{i}=vec(F==u(i)); F_rep(i)=F_idx{i}(1,1); end; clear i vec l u;

%%%%%%%%%%%%%%%%%%%%%%%%below is fixed %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%counterfactual values
CEn=zeros(N,1); Cpr=zeros(N,1); Cq=zeros(N,1); CCS=zeros(N,1);
Cpmc=zeros(N,1); Cmarkup=zeros(N,1); NEn=zeros(N,1);
CX=[En Ex NR NXFD];
%%%%%%%%%%%% Restrict the data to year 2010 %%%%%%%%%%%%%%%%%%
%%Run estimation for every quarter in the year %%%
for iQ=1:4
S=find(Y==2010&Q==iQ); %index for counterfactual parts
%%% Only care about the experiment period S.
%Vectors starting with S have dimension Sx1 %%
Sc_cjt=c_cjt(S); SM=M(S);
STC=TC(S); STF=TF(S); SEn=En(S); SN=size(STC,1);
SX=CX(S,:); KSX=size(SX,2); Slocal=local(S,:); 
beta1=betaX(1); KbetaX=size(betaX,1);
eX2beta2=exp(SX(:,2:KSX)*betaX(2:KbetaX));
%%% Get a representative index in each group in TJ
u=unique(STF,'stable'); L=length(u); STF_idx=cell(L,1); vec=1:SN; STF_rep=zeros(1,L);
for i=1:L, STF_idx{i}=vec(STF==u(i)); STF_rep(i)=STF_idx{i}(1,1); end;
clear i vec l u;
%%% Get a representative index in each group in TJ
u=unique(STC,'stable'); L=length(u); STC_idx=cell(L,1); vec=1:SN; STC_rep=zeros(1,L);
for i=1:L, STC_idx{i}=vec(STC==u(i)); STC_rep(i)=STC_idx{i}(1,1); end; clear i vec l u;
%%Shorten price and cost to TJx1, instead of CTJx1
P_jt=SEn(STF_rep); Sc_jt=Sc_cjt(STF_rep);
%%Define the equation to solve for equilibrium price
fun_price_importban=@(P_jt)solveprice_importban(beta1,P_jt,eX2beta2,Sc_jt,c,SM,STC_idx,STF_idx,STF_rep,SN,mM);
opfsolve=optimoptions('fsolve','Algorithm','levenberg-marquardt','Display','iter',...
    'MaxIterations',400,'MaxFunctionEvaluations',1e15,...
    'StepTolerance',1e-10,'FunctionTolerance',1e-10,'OptimalityTolerance',1e-6); 
%'Algorithm','levenberg-marquardt'
%%Calculate eq price for trash tax
[SP,feval,exitflag,output,jacob]=fsolve(fun_price_importban,P_jt,opfsolve);
CEn(intersect(TF_rep,S))=SP; %P_jt in period S
CEn(S)=my_accumarray(STF_idx,CEn(S));%span. same as egen x=mean(.) or total
SX=[CEn(S) SX(:,2:KSX)];
[Cpr(S),CCS(S)]=pr_fun(betaX,SX,STC_idx); %counterfactual market share
%% calculate new price cost margin
dsdp(S)=betaX(1)*Cpr(S).*(1-Cpr(S)); %dsj/dpj
nume(S)=SM./mean(M).*Cpr(S);
nume(S)=my_accumarray(STF_idx,nume(S));
deno(S)=my_accumarray(STF_idx,nume(S).*dsdp(S));
Cpmc(S)=-nume(S)./deno(S); Cmarkup(S)=Cpmc(S)./(CEn(S)-Cpmc(S));
end
%out of the loop S
Cq=Cpr.*M; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=CEn; %new eq price. Nrealprice=CEn*tax...
mysample.Cpmc=Cpmc; mysample.Cmarkup=Cmarkup;
writetable(mysample,strcat('output\nimby_counterfactual_importBan_with_eqprice_',counter(coun),'.xlsx'),'Sheet',1);

%%%%%%%%%%%%% IMPORT TAX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(strcat('output\result_main_estimator'));
M=sj_q; mM=mean(M);
cnty=mysample.desco==counter(coun);
%%%%%%%%%%% Counterfactuals %%%%%%%%%%%%%%%%%%%%%%%%%%%
%Claim variables and define counterfactual for fuel tax and import ban;
%For trash tax and import tax, modify the function that calculate eq price
%keep local observations for mysample. Define obs to remove
%i.e. remove nonlocalf and not 2010
remove=mysample.year~=2010;% zeros(N,1);
mysample(remove==1,:)=[]; %only keep local flow sample
c_cjt(remove==1,:)=[];
M(remove==1,:)=[];
cnty(remove==1,:)=[];
NXFD=XFD(remove==0,:); NR=R(remove==0,:);
%redefine variables from mysample
En=mysample(:,'realprice'); En=table2array(En);
Ex=mysample(:,{'dridis','localf'}); Ex=table2array(Ex); KEx=size(Ex,2);
T=mysample(:,'T');T=table2array(T);
Y=mysample(:,'year');Y=table2array(Y);
Q=mysample(:,'quarter');Q=table2array(Q);
TC=mysample(:, 'TC'); TC=table2array(TC); TCmax=max(TC);
C=mysample(:,'C'); C=table2array(C); 
F=mysample(:,'F'); F=table2array(F);
TF=mysample(:,'TF'); TF=table2array(TF);
q=mysample(:,'q'); q=table2array(q); 
local=mysample(:,'localf'); local=table2array(local);
q1=q>0; q0=q==0; lq=log(q); lq(q0)=0;
N=size(q,1); %number of observations
%%%% Get obs indices of each group in TC. Each cell is a tc group
uTC=unique(TC,'stable'); L=length(uTC); TC_idx=cell(L,1); vec=1:N; TC_rep=zeros(1,L);
TC_l=NaN(N,1);
for i=1:L
    %all obs indices of group tc=u(i). each cell is for 1 group
    TC_idx{i}=vec(TC==uTC(i)); 
    %representative indices of all groups, length(v)=length(u)
    %first representative index of each group tc=u(i)
    TC_rep(i)=TC_idx{i}(1,1); 
    TC_l(TC_idx{i})=size(TC_idx{i},2);
end
clear i vec l;
TCmj_idx=cell(N,1); %obs indices of k in TC_{j} but exclude facility j
for i=1:N, TCmj_idx{i}=setdiff(TC_idx{uTC==TC(i)},i); end
%%% Get a representative index in each group in T
u=unique(T,'stable'); L=length(u); T_idx=cell(L,1); vec=1:N; T_rep=zeros(1,L);
for i=1:L, T_idx{i}=vec(T==u(i)); T_rep(i)=T_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in C
u=unique(C,'stable'); L=length(u); C_idx=cell(L,1); vec=1:N; C_rep=zeros(1,L);
for i=1:L, C_idx{i}=vec(C==u(i)); C_rep(i)=C_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in TF
u=unique(TF,'stable'); L=length(u); TF_idx=cell(L,1); vec=1:N; TF_rep=zeros(1,L);
for i=1:L, TF_idx{i}=vec(TF==u(i)); TF_rep(i)=TF_idx{i}(1,1); end; clear i vec l u;
%%% Get a representative index in each group in F
u=unique(F,'stable'); L=length(u); F_idx=cell(L,1); vec=1:N; F_rep=zeros(1,L);
for i=1:L, F_idx{i}=vec(F==u(i)); F_rep(i)=F_idx{i}(1,1); end; clear i vec l u;

%%%%%%%%%%%%%%%%%%%%%%%%below is fixed %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CEn=zeros(N,1); Cpr=zeros(N,1); Cq=zeros(N,1); CCS=zeros(N,1);
Cpmc=zeros(N,1); Cmarkup=zeros(N,1); NEn=zeros(N,1);
CX=[En Ex NR NXFD];
%%%%%%%%%%%% Restrict the data to year 2010 %%%%%%%%%%%%%%%%%%
%%Run estimation for every quarter in the year %%%
for iQ=1:4
S=find(Y==2010&Q==iQ); %define counterfactual period
Scnty=cnty(S);
%%% Only care about the experiment period S.
%Vectors starting with S have dimension Sx1 %%
Sc_cjt=c_cjt(S); SM=M(S)./1; 
STC=TC(S); STF=TF(S); SEn=En(S); SN=size(STC,1);
SX=CX(S,:); KSX=size(SX,2); Slocal=local(S,:);
beta1=betaX(1); KbetaX=size(betaX,1);
eX2beta2=exp(SX(:,2:KSX)*betaX(2:KbetaX));
%%% Get a representative index in each group in TJ
u=unique(STF,'stable'); L=length(u); STF_idx=cell(L,1); vec=1:SN; STF_rep=zeros(1,L);
for i=1:L, STF_idx{i}=vec(STF==u(i)); STF_rep(i)=STF_idx{i}(1,1); end;
clear i vec l u;
%%% Get a representative index in each group in TJ
u=unique(STC,'stable'); L=length(u); STC_idx=cell(L,1); vec=1:SN; STC_rep=zeros(1,L);
for i=1:L, STC_idx{i}=vec(STC==u(i)); STC_rep(i)=STC_idx{i}(1,1); end; clear i vec l u;
%%Shorten price and cost to TJx1, instead of CTJx1
P_jt=SEn(STF_rep); Sc_jt=Sc_cjt(STF_rep);
%old price/observed price is guessing point
%%Define the equation to solve for equilibrium price
fun_price_importtax=@(P_jt)solveprice_importtax_cnty(beta1,P_jt,eX2beta2,Sc_jt,c,SM,Slocal,Scnty,STC_idx,STF_idx,STF_rep,SN,mM);
opfsolve=optimoptions('fsolve','Algorithm','levenberg-marquardt','Display','iter',...
    'MaxIterations',400,'MaxFunctionEvaluations',1e15,...
    'StepTolerance',1e-12,'FunctionTolerance',1e-12,'OptimalityTolerance',1e-5); 
%%Calculate eq price for trash tax
[SP,feval,exitflag,output,jacob]=fsolve(fun_price_importtax,P_jt,opfsolve);
%span P_jt to P_cjt
CEn(intersect(TF_rep,S))=SP; %P_jt;
CEn(S)=my_accumarray(STF_idx,CEn(S));
NEn(S)=CEn(S);
NEn(S & local(S)==0&Scnty==1)=CEn(S & local(S)==0&Scnty==1)*1.10; % price includes tax
SX=[NEn(S) SX(:,2:KSX)];
[Cpr(S),CCS(S)]=pr_fun(betaX,SX,STC_idx); %counterfactual market share
%% calculate new price cost margin
dsdp(S)=betaX(1)*Cpr(S).*(1-Cpr(S)); %dsj/dpj
nume(S)=SM./mean(M).*Cpr(S);
nume(S)=my_accumarray(STF_idx,nume(S));
deno(S)=my_accumarray(STF_idx,nume(S).*dsdp(S));
Cpmc(S)=-nume(S)./deno(S); Cmarkup(S)=Cpmc(S)./(CEn(S)-Cpmc(S));
end
%out of the loop S
Cq=Cpr.*M; %counterfactual trash amount in a facility
%%write into csv to input STATA for graphical analysis
mysample.Cpr=Cpr;
mysample.Cq=Cq;
mysample.CCS=CCS;
mysample.Crealprice=CEn; 
mysample.Cpmc=Cpmc; mysample.Cmarkup=Cmarkup;
writetable(mysample,strcat('output\nimby_counterfactual_importTax_with_eqprice_',counter(coun),'.xlsx'),'Sheet',1);

end %for loop in county for counterfactual

%%%%%%%%% END HERE. %%%%%%%%%%%%%%%%%%%%%%%%%
diary off